The present invention relates to processing of seismic data representative of subsurface features in the earth and, more particularly, to improved methods of processing seismic data to remove unwanted noise from meaningful reflection signals.
Seismic surveys are one of the most important techniques for discovering the presence of oil and gas deposits. If the data is properly processed and interpreted, a seismic survey can give geologists a picture of subsurface geological features, so that they may better identify those features capable of holding oil and gas. Drilling is extremely expensive, and ever more so as easily tapped reservoirs are exhausted and new reservoirs are harder to reach. Having an accurate picture of an area""s subsurface features can increase the odds of hitting an economically recoverable reserve and decrease the odds of wasting money and effort on a nonproductive well.
The principle behind seismology is deceptively simple. As seismic waves travel through the earth, portions of that energy are reflected back to the surface as the energy waves traverse different geological layers. Those seismic echoes or reflections give valuable information about the depth and arrangement of the formations, some of which hopefully contain oil or gas deposits.
A seismic survey is conducted by deploying an array of energy sources and an array of sensors or receivers in an area of interest. Typically, dynamite charges are used as sources for land surveys, and air guns are used for marine surveys. The sources are discharged in a predetermined sequence, sending seismic energy waves into the earth. The reflections from those energy waves or xe2x80x9csignalsxe2x80x9d then are detected by the array of sensors. Each sensor records the amplitude of incoming signals over time at that particular location. Since the physical location of the sources and receivers is known, the time it takes for a reflection signal to travel from a source to a sensor is directly related to the depth of the formation that caused the reflection. Thus, the amplitude data from the array of sensors can be analyzed to determine the size and location of potential deposits.
This analysis typically starts by organizing the data from the array of sensors into common geometry gathers. That is, data from a number of sensors that share a common geometry are analyzed together. A gather will provide information about a particular spot or profile in the area being surveyed. Ultimately, the data will be organized into many different gathers and processed before the analysis is completed and the entire survey area mapped.
The types of gathers typically used include: common midpoint, where the sensors and their respective sources share a common midpoint; common source, where all the sensors share a common source; common offset, where all the sensors and their respective sources have the same separation or xe2x80x9coffsetxe2x80x9d; and common receiver, where a number of sources share a common receiver. Common midpoint gathers are the most common gather today because they allow the measurement of a single point on a reflective subsurface feature from multiple source-receiver pairs, thus increasing the accuracy of the depth calculated for that feature.
The data in a gather is typically recorded or first assembled in the offset-time domain. That is, the amplitude data recorded at each of the receivers in the gather is assembled or displayed together as a function of offset, i.e., the distance of the receiver from a reference point, and as a function of time. In theory, when the gathered data is displayed in the offset-time domain, or xe2x80x9cX-Txe2x80x9d domain, the amplitude peaks corresponding to reflection signals detected at the gathered sensors should align into patterns that are referred to as kinematic travel time trajectories. It is from those trajectories that one ultimately may determine an estimate of the depths at which formations exist.
A number of factors, however, make the practice of seismology and, especially, the interpretation of seismic data much more complicated than its basic principles. First, the reflected signals that indicate the presence of geological strata typically are mixed with a variety of noise.
The most meaningful signals are the so-called primary reflection signals, those signals that travel down to the reflective surface and then back up to a receiver. When a source is discharged, however, a portion of the signal travels directly to receivers without reflecting off of any subsurface features. In addition, a signal may bounce off of a subsurface feature, bounce off the surface, and then bounce off the same or another subsurface feature, one or more times, creating so-called multiple reflection signals. Other portions of the signal turn into noise as part of ground roll, refractions, and unresolvable scattered events. Some noise, both random and coherent, is generated by natural and man-made events outside the control of the survey.
All of this noise is occurring simultaneously with the reflection signals that indicate subsurface features. Thus, the noise and reflection signals tend to overlap when the survey data is displayed in X-T space. The overlap can mask primary reflection signals and make it difficult or impossible to identify patterns in the display upon which inferences about subsurface geology may be drawn. Accordingly, various mathematical methods have been developed to process seismic data in such a way that noise is separated from primary reflection signals.
Many such methods seek to achieve a separation of signal and noise by transforming the data from the X-T domain to other domains. In other domains, such as the frequency-wavenumber (F-K) domain or the time-slowness (Tau-P), there is less overlap between the signal and noise data. Once the data is transformed, various mathematical filters are applied to the transformed data to eliminate as much of the noise as possible and, thereby, to enhance the primary reflection signals. The data then is inverse transformed back into the offset-time domain for interpretation or further processing.
For example, so-called Radon filters are commonly used to attenuate or remove multiple reflection signals. Such methods typically first process a data gather to compensate for the increase in travel time as sensors are further removed from the source (normal moveout or xe2x80x9cNMOxe2x80x9d correction). After NMO correction, the data is transformed along kinematic travel time trajectories in offset-time space to corresponding trajectories in time-slowness space, where slowness p is defined as reciprocal velocity (or p=1/xcexd). Primary reflection signal trajectories will transform into different regions of time-slowness space than the multiple reflection signals. Thus, the primary reflection signals are isolated from multiple reflections and may be muted, leaving only the signals for multiple reflections. The data is then transformed back into offset-time space, and the data which has not been muted is then subtracted from the original data gather. Thus, the multiple reflection signals are removed from the data gather, leaving the primary reflection signals more readily apparent and easier to interpret.
It will be appreciated, however, that in such prior art Radon filters, the multiples recorded by receivers closer to the gather reference point are not as effectively separated from the primary reflections. Thus, those multiples are not subtracted from the original data gather and, when the data ultimately is interpreted, remain present as noise that can mask primary reflection signals.
Radon filters have been developed for use in connection with common source, common receiver, common offset, and common midpoint gathers. They include those based on linear slant-stack, parabolic, and hyperbolic kinematic travel time trajectories. The general case forward transformation equation used in Radon filtration processes, R(p, xcfx84)[d(x, t)], is set forth below:       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢          xe2x80x83        ⁢                  ⅆ        x            ⁢                        ∫                      -            ∞                    ∞                ⁢                  xe2x80x83                ⁢                              ⅆ            t                    ⁢                      xe2x80x83                    ⁢                      ⅆ                          (                              x                ,                t                            )                                ⁢                      δ            ⁡                          [                              f                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                          
(forward transformation)
where
u(p, xcfx84)=transform coefficient at slowness p and zero-offset time
xcfx84d(x, t)=measured seismogram at offset x and two-way time t
p=slowness
t=two-way travel time
xcfx84=two-way intercept time at p=0
x=offset
xcex4=Dirac delta function
ƒ(t, x, xcfx84, p)=forward transform function
The forward transform function for linear slant stack kinematic travel time trajectories is as follows:
ƒ(t, x, xcfx84, p)=txe2x88x92xcfx84xe2x88x92px
where
xcex4[ƒ(t, x, xcfx84, p)]=xcex4(txe2x88x92xcfx84xe2x88x92px)
xe2x80x83=1, when t=xcfx84+px, and
xe2x80x83=0, elsewhere.
Thus, the forward linear slant stack Radon transformation equation becomes       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢              xe2x80x83            ⁢              ⅆ                  (                      x            ,                          τ              +              px                                )                    
The forward transform function for parabolic kinematic trajectories is as follows:
ƒ(t, x, xcfx84, p)=txe2x88x92xcfx84xe2x88x92px2
where
xcex4[ƒ(t, x, xcfx84, p)]=xcex4(txe2x88x92xcfx84xe2x88x92px2)
xe2x80x83=1, when t=xcfx84+px2, and
xe2x80x83=0, elsewhere.
Thus, the forward parabolic Radon transformation equation becomes       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢              xe2x80x83            ⁢              ⅆ                  (                      x            ,                          τ              +                              px                2                                              )                    
The forward transform function for hyperbolic kinematic travel time trajectories is as follows:
ƒ(t, x, xcfx84, p)=txe2x88x92{square root over (xcfx842+p2x2)}
where   "AutoLeftMatch"                                          δ            ⁡                          [                              f                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                =                      δ            ⁡                          (                              t                -                                                                            τ                      2                                        +                                                                  p                        2                                            ⁢                                              x                        2                                                                                                        )                                                                                =            1                    ,                                    when              ⁢                              xe2x80x83                            ⁢              t                        =                                                            τ                  2                                +                                                      p                    2                                    ⁢                                      x                    2                                                                                ,          and                                                          =            0                    ,                      elsewhere            .                              
Thus, the forward hyperbolic Radon transformation equation becomes       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢              xe2x80x83            ⁢              ⅆ                  (                      x            ,                                                            τ                  2                                +                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
A general case inverse Radon transformation equation is set forth below:       d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢          xe2x80x83        ⁢                  ⅆ        p            ⁢                        ∫                      -            ∞                    ∞                ⁢                  xe2x80x83                ⁢                              ⅆ            τ                    ⁢                      xe2x80x83                    ⁢                                    ρ              ⁡                              (                τ                )                                      *                    ⁢                      u            ⁡                          (                              p                ,                τ                            )                                ⁢                      xe2x80x83                    ⁢                      δ            ⁡                          [                              g                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                          
(inverse transformation)
where
g(t, x, xcfx84, p)=inverse transform function, and
xcfx81(xcfx84)*=convolution of rho filter.
The inverse transform function for linear slant stack kinematic trajectories is as follows:
g(t, x, xcfx84, p)=xcfx84xe2x88x92t+px
Thus, when xcfx84=txe2x88x92px, the inverse linear slant stack Radon transformation equation becomes       d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢          xe2x80x83        ⁢                  ⅆ        p            ⁢              xe2x80x83            ⁢                        ρ          ⁡                      (            τ            )                          *            ⁢              u        ⁡                  (                      p            ,                          t              -              px                                )                    
The inverse transform function for parabolic trajectories is as follows:
g(t, x, xcfx84, p)=xcfx84xe2x88x92t+px2
Thus, when xcfx84=txe2x88x92px2, the inverse linear slant stack Radon transformation equation becomes       d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        p            ⁢              xe2x80x83            ⁢      ρ      ⁢              xe2x80x83            ⁢                        (          τ          )                *            ⁢              u        ⁡                  (                      p            ,                          t              -                              px                2                                              )                    
The inverse transform function for hyperbolic trajectories is as follows:
g(t, x, xcfx84, p)=xcfx84xe2x88x92{square root over (t2xe2x88x92p2x2)}
Thus, when xcfx84={square root over (t2xe2x88x92p2x2)}, the inverse hyperbolic Radon transformation equation becomes       d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        p            ⁢              xe2x80x83            ⁢      ρ      ⁢              xe2x80x83            ⁢                        (          τ          )                *            ⁢              u        ⁡                  (                      p            ,                                                            t                  2                                -                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
To date, Radon transformations based on linear slant stack and parabolic trajectories have been more commonly used. In common midpoint gathers, however, primary reflection signals will generally have hyperbolic kinematic trajectories, thus, it would be preferable to use hyperbolic Radon transforms. The computational intensity of those equations, however, typically has precluded the use of fine sampling variables in performing the transformation.
In particular it will be noted that the transform function for hyperbolic trajectories contains a square root whereas those for linear slant stack and parabolic transform functions do not. Thus, as compared to linear slant stack and parabolic Radon transformations, the sampling variables used in performing hyperbolic Radon transformations are relatively coarse to compensate for the additional computations needed to handle square root computations. Increasing the magnitude of the sampling variables reduces the number of computations needed to process data, but the relative coarseness of the sampling variables reduces the accuracy of the process.
Forward and inverse hyperbolic Radon transforms also are not exact inverses of each other. Accordingly, an additional least-mean-squares (xe2x80x9cLMSxe2x80x9d) energy minimization calculation is often used. Also, the data is typically processed to compensate for the diminution in amplitude as a signal travels through the earth (xe2x80x9camplitude balancingxe2x80x9d), to NMO correct, or for various other reasons that require additional processing. Those calculations in general and, especially LMS analyses, are complex and require substantial computing time, even on advanced computer systems.
Moreover, a typical seismic survey may involve hundreds of sources and receivers, thus generating tremendous quantities of raw data. The data may be subjected to thousands of different data gathers. Once the gathers have been filtered, the gathers must be stacked, or compiled with other gathers, and subjected to further processing in order to make inferences about subsurface formations. Thus, and especially so when conventional hyperbolic Radon transformations and LMS calculations are used, the cost of processing seismic data can be extremely high. Moreover, because of the computational intensity of prior art processes and the attendant cost, it has been necessary to avoid the use of hyperbolic Radon transformations or to use relatively coarse sampling variables, thereby saving computations and cost. The computational savings, however, come at the price of reduced accuracy.
The velocity at which primary reflection signals are traveling, what is referred to as the stacking velocity, is used in various analytical methods that are applied to seismic data. For example, it is used in determining the depth and lithology or sediment type of geological formations in the survey area. It also is used various seismic attribute analyses. A stacking velocity function may be determined by what is referred to as a semblance analysis. Which such methods have provided reasonable estimations of the stacking velocity function for a data set, given the many applications based thereon, it is desirable to define the stacking velocity function as precisely as possible.
An object of this invention, therefore, is to provide a method for processing seismic data that more effectively removes unwanted noise from meaningful reflection signals.
It also is an object to provide such methods based on hyperbolic Radon transformations.
Another object of this invention is to provide methods for removal of noise from seismic data that are comparatively simple and require relatively less computing time and resources.
Yet another object is to provide methods for more accurately defining the stacking velocity function for a data set.
It is a further object of this invention to provide such methods wherein all of the above-mentioned advantages are realized.
Those and other objects and advantages of the invention will be apparent to those skilled in the art upon reading the following detailed description and upon reference to the drawings.
The subject invention provides for methods of processing seismic data to remove unwanted noise from meaningful reflection signals indicative of subsurface formations. The methods comprise the steps of obtaining field records of seismic data detected at a number of seismic receivers in an area of interest. The seismic data comprises signal amplitude data recorded over time by each of the receivers and contains both primary reflection signals and unwanted noise events. The seismic data is assembled into common geometry gathers in an offset-time domain.
An offset weighting factor xn is then applied to the amplitude data, wherein 0 less than n less than 1, and the offset weighted data is transformed from the offset-time domain to the time-slowness domain using a Radon transformation. Such Radon transformations include the following continuous transform equation, and discrete versions thereof that approximate the continuous transform equation:       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢          xe2x80x83        ⁢                  ⅆ        x            ⁢                        ∫                      -            ∞                    ∞                ⁢                  xe2x80x83                ⁢                              ⅆ            t                    ⁢                      xe2x80x83                    ⁢                      ⅆ                          (                              x                ,                t                            )                                ⁢                      x            n                    ⁢                      xe2x80x83                    ⁢                      δ            ⁡                          [                              f                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                          
Radon transformations based on linear slant stack, parabolic, and other non-hyperbolic kinematic travel time trajectories may be used, but those based on hyperbolic Radon kinematic trajectories are preferred. Suitable hyperbolic Radon transformations include the following continuous transform equation, and discrete versions thereof that approximate the continuous transform equation:       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢          xe2x80x83        ⁢                  ⅆ        x            ⁢              xe2x80x83            ⁢              x        n            ⁢              ⅆ                  (                      x            ,                                                            τ                  2                                +                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
Preferably, only the datum having a slowness pd that is greater that a predetermined minimum, pmin, and less than a predetermined maximum, pmax, is transformed.
A corrective filter is then applied to enhance the primary reflection signal content of the data and to eliminate unwanted noise events. Preferably, the corrective filter is determined by performing a semblance analysis on the amplitude data to generate a semblance plot. A velocity analysis is then performed on the semblance plot, from which the corrective filter is defined. Preferably, the corrective filter is a time variant velocity filter.
After filtering, the enhanced signal content is inverse transformed from the time-slowness domain back to the offset-time domain using an inverse Radon transformation, and an inverse of the offset weighting factor pn is applied to the inverse transformed data, wherein 0 less than n less than 1. Such Radon transformations include the following continuous inverse transform equation, and discrete versions thereof that approximate the continuous inverse transform equation:       d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢          xe2x80x83        ⁢                  ⅆ        p            ⁢                        ∫                      -            ∞                    ∞                ⁢                  xe2x80x83                ⁢                              ⅆ            τ                    ⁢                      xe2x80x83                    ⁢                      p            n                    ⁢                      xe2x80x83                    ⁢                                    ρ              ⁡                              (                τ                )                                      *                    ⁢                      u            ⁡                          (                              p                ,                τ                            )                                ⁢                      xe2x80x83                    ⁢                      δ            ⁡                          [                              g                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                          
Suitable hyperbolic inverse Radon transformations include the following continuous inverse transform equation, and discrete versions thereof that approximate the continuous inverse transform equation:       d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        p            ⁢              xe2x80x83            ⁢              p        n            ⁢              xe2x80x83            ⁢      ρ      ⁢              xe2x80x83            ⁢                        (          τ          )                *            ⁢              u        ⁡                  (                      p            ,                                                            t                  2                                -                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
The signal amplitude data for the enhanced signal content is thereby restored. The processed and filtered data may then be subject to further processing by which inferences about the subsurface geology of the survey area may be made.
It will be appreciated that primarily because they utilize an offset weighting factor and its inverse, and avoid more complex mathematical operations required by prior art processes, the novel methods require less total computation time and resources. Thus, even though hyperbolic Radon transformations heretofore have generally been avoided because of their greater complexity relative to linear slant stack and parabolic Radon transformations, the novel processes are able to more effectively utilize hyperbolic Radon transformations and take advantage of the greater accuracy they can provide, especially when applied to common midpoint geometry gathers.
At the same time, because they are computationally less intensive, the sampling variables used in the transformations and semblance analyses may be set much finer than are typically used in prior art processes. xcex94t and xcex94xcfx84 values may be set as low as the sampling rates at which the data was recorded, and xcex94p values as low as about 0.5 xcexcsec/m may be used. The use of finer sampling variables also allows for more accurate processing of data.
The subject invention also provides for methods of filtering seismic data to remove unwanted noise from meaningful reflection signals indicative of subsurface formations that comprise obtaining field records of seismic data detected at a number of seismic receivers in an area of interest. The seismic data comprises signal amplitude data recorded over time by each receiver and contains both primary reflection signals and unwanted noise events. Each of the amplitude datum has an associated slowness pd.
The data is assembled into common geometry gathers in an offset-time domain. The amplitude data where pd is greater that a predetermined minimum, pmin, and less than a predetermined maximum, pmax, is then transformed from the offset-time domain to the time-slowness domain using a Radon transformation.
A corrective filter is then applied to the transformed data to enhance the primary reflection signal content of the data and to eliminate unwanted noise events, and the enhanced signal content is inverse transformed from the time-slowness domain back to the offset-time domain using an inverse Radon transformation. The signal amplitude data for the enhanced signal content is thereby restored.
Other methods of the subject invention comprise assembling seismic amplitude data into common geometry gathers in an offset-time domain and transforming the amplitude data from the offset-time domain to the time-slowness domain using a Radon transformation.
A semblance analysis is also performed on the amplitude data to generate a semblance plot. A velocity analysis is then performed on the semblance plot to define a stacking velocity function and a time variant corrective filter to enhance the primary reflection signal content of the transformed data and to eliminate unwanted noise events.
The corrective filter is then applied to the transformed data. The filtered date is inverse transformed from the time-slowness domain back to the offset-time domain using an inverse Radon transformation, thereby restoring the signal amplitude data for the filtered data.
Other methods comprise designing a Tau-P transform for each or the X-T common geometry gathers into which the seismic data has been assembled to determine the number of slownesses p and the xcex94p value. A given choice of numbers of p""s are then applied as a family of hyperbolas over the CMP gather""s X-T space. Amplitudes over the family of hyperbolas defined by each p value are summed as a constant velocity hyperbolic time trajectory, and the amplitude data is transformed from the offset-time domain to the time-slowness domain using a Radon transformation.
A semblance analysis is then performed on the amplitude data using the xcex94p value to generate a semblance plot. A velocity analysis is performed on the semblance plot to define a stacking velocity function and a time variant velocity filter to enhance the primary reflection signal content of the transformed data and to eliminate unwanted noise events, and the filter is applied to the transformed data. The filtered data is inverse transformed from the time-slowness domain back to the offset-time domain using an inverse Radon transformation, thereby restoring the signal amplitude data for the filtered data.
The subject invention also provides for methods for determining a stacking velocity function. The methods comprise the steps of processing and filtering the data as described above. A semblance analysis is performed on the restored data to generate a refined semblance plot, and a velocity analysis is performed on the refined semblance plot to define a refined stacking velocity function. Preferably, an offset weighting factor xn is first applied to the restored amplitude data, wherein 0 less than n less than 1.
It will be appreciated that, because the semblance and velocity analysis is performed on data that has been more accurately processed and filtered, that the resulting stacking velocity function is more accurate. This enhances the accuracy of the many different seismic processing methods that utilize stacking velocity functions. Ultimately, the increased accuracy and efficiency of the novel processes enhances the accuracy of surveying underground geological features and, therefore, the likelihood of accurately locating the presence of oil and gas deposits.